Charmonium mass in hot and dense hadronic matter 
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We investigate mass shifts of charmonia driven by change of the gluon condensate below but near 
transition temperatures at finite baryonic chemical potential. Extending previous prescription on the 
relation between gluon condensates and thermodynamic quantities, we model the gluon condensates 
of hadronic matter at finite temperature and baryonic chemical potential such that the scalar gluon 
condensate fits with the latest lattice QCD data. By making use of the QCD sum rule and the second 
order Stark effect, we find that the smoother transition in the full QCD can lead to moderate mass 
shifts of charmonia even below the transition temperature. We also find larger mass shift at fixed 
temperature as chemical potential increases. Existing data on charmonium-charmonium ratio is 
found to be consistent with the statistical hadronization scenario including the obtained mass shift. 

PACS numbers: 14.40.Pq,11.55.Hx,24.85.+p 



I. INTRODUCTION 

Properties of heavy quarkonia in medium have been 
extensively studied since it was pointed out that sup- 
pression of J I ip by Debye screening could be a signature 
of creation of the deconfined matter in relativistic heavy 
ion collisions [l|. It should be noted that a mass shift of 
a charmonium state in hot hadronic environment could 
be a precursor phenomenon of the transition caused by 
a decrease in the string tension @ . These early expecta- 
tions are based on an intuitive picture on a quarkonium, 
a heavy quark and its antiquark bound by a confining 
potential, which successfully describes the properties in 
vacuum [3j . While the lattice QCD provides a first prin- 
ciple approach to the problem, it still lacks the necessary 
resolution needed for discriminating possible changes in 
the quarkonium spectral function at finite temperature, 
especially near the critical temperature where an abrupt 
change could take place. The maximum entropy method 
for this problem can at best only tell us about the (non- 
existence of the lowest peak in the spectral function 0- 
8]. Therefore, to assess the medium modification one 
needs a complementary framework such as the poten- 
tial model which utilizes a quark- antiquark potential ex- 
tracted from lattice calculation In the meantime, we 
have proposed an approach utilizing the QCD sum rule 
and the second order Stark effect which allows us to re- 
late the temperature dependent gluon condensates as the 
primary inputs from lattice QCD to the spectral changes 
of heavy quarkonia [Iol-[l4| . 

So far our studies relied on the gluon condensates ex- 
tracted from the trace anomaly of pure SU(3) case (l5j . 
In pure gauge theory with N c > 3, there is a first or- 
der deconfinement transition at T = T c thus the trace 



'Electronic address: kmorita@yukawa.kyoto-u.ac.jp 
^Electronic address: suhoung@yonsei.ac.kr 



anomaly shows abrupt change across T c [161 ] , which leads 
to the similar behavior of the spectral property of the 
quarkonia. 

In order to compare results with experimental data, we 
need more realistic estimates of the condensates based on 
full QCD lattice calculations including dynamical light 
quarks. Recently the calculation of the equation of state 
has been carried out with physical quark masses at van- 
ishing chemical potential [17j . Due to the crossover na- 
ture of the transition [l8j , the critical temperature is no 
longer a well-defined quantity. The pseudocritical tem- 
perature defined by a peak or an inflection point depends 
on observables [19| and is found to range from 147 MeV 
(chiral susceptibility) to 165 MeV (strange quark number 
susceptibility). 

As for the scalar gluon condensate, which is the glu- 
onic part of the trace anomaly, it is expected to have the 
same bulk property but with smoother change near the 
pseudocritical temperature T pc . In fact, the magnitude 
for the change of the scalar condensate in a pure gauge 
theory at the transition region has been found to be al- 
most the same as in the full QCD case [ll|, H3| when the 
temperatures are normalized by T c and T pc , despite the 
difference in the (pseudo)critical temperatures; T c ~ 265 
MeV in the pure SU(3) [l![l]t and T pc ~ 190-200 MeV 
in full (2+1) flavor QCD with m ff ~ 220 MeV 0. The 
crossover nature of the transition has led to the smoother 
tem per ature dependence of the gluon condensate near 
Tp C [2lJ. For the twist- 2 gluon condensate, which is not 
a dominant but non-negligible contribution to the sum 
rule [HI H2] , no lattice data at physical quark masses is 
available yet. 

In this paper, we first extract the gluon condensates in 
full QCD by making use of the resonance gas model. This 
prescription enables us to study finite baryonic chemical 
potential case also. Using these gluon condensates, we 
investigate the change of spectral properties of charmo- 
nia in hot hadronic matter at various values of temper- 
ature and baryonic chemical potential along the freeze- 
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out line which is deduced from statistical model analyses 
[23l [2^ and is expected to be close to the hadroniza- 
tion points. Although the charmonium production mech- 
anism in rclativistic heavy ion collisions has not been 
understood well, some experimental data seem to indi- 
cate that the statistical production of charmonia could 
be possible [25|. Despite the complexity of the colli- 
sion processes, this scenario considerably simplifies the 
charmonium-charmonium particle number ratio. We ex- 
amine possible influences of the spectral modification on 
this quantity. 

In the next section, we present a resonance gas model 
for the gluon condensates of hot and dense hadronic mat- 
ter. In Sec Hm we report results of spectral changes of 
charmonia. The experimental implication will be dis- 
cussed in Sec.|lVl Section [V] is devoted to a summary. 



II. RESONANCE GAS MODEL FOR THE 
GLUON CONDENSATES 

We start with two quantities Mq and M2 characterizing 
the temperature dependence of the thermal expectation 
value of gluonic operators 



(1) 



u a u p - -g afs ) M 2 (T). (2) 



Here, the symbol ST denotes the traceless and symmet- 
ric part of the operator. We assign to Mq(T) only the 
temperature dependent part of the expectation value; it 
has in general a temperature independent part which is 
nothing but the gluon condensate in the vacuum. 

The above equations immediately relate Mq and M2 to 
the thermodynamic quantities via the energy-momentum 
tensor in thermal equilibrium; namely, Mq(T) = e — 3p 
and M 2 (T) = e+p with e and p being the energy density 
and the pressure respectively in the case of pure gluonic 
system jl 31 ] - In the presence of fermions, however, it is 
not straightforward to relate Mq and M2 to the thermo- 
dynamic quantities since the energy-momentum tensor 
has fermionic part. The trace anomaly receives contribu- 
tions from massive fermions as 
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^2m l {q i q i ). 



(3) 



Mq becomes e — 3p only when the second term is negli- 
gible, otherwise, the fermionic part has to be explicitly 
subtracted out from e — 3p before identifying it to Mq. In 
lattice QCD calculations, contributions from each term 
in Eq. © have been estimated [13, HJ 1 . Similarly, M 2 



1 In Ref. jl7| |. however, somewhat different scheme is used to cal- 
culate the trace anomaly. 



can be related to the off-diagonal part of the energy- 
momentum tensor after the fermionic part is subtracted 
out. Such data are not yet available. Therefore we need 
a scheme to subtract the fermionic contribution from the 
total of the energy-momentum tensor. 

Going back to the original definition given in Eq. ©, 
one can relate the nucleon expectation values to Mq and 
Mi within the linear density approximation 26] 



M n - m - 
M „.m. 



pm%, 
pA G m N , 



(4) 
(5) 



where p, m° N , Aq and mjv are the density of nucleus, the 
nucleon mass in the chiral limit, the second moment of 
gluon distribution function of the nucleon, and the nu- 
cleon mass, respectively. One sees in Eq. ^ that the 
chiral limit is taken for the nucleon mass, which cor- 
responds to removing the fermionic term in the trace 
anomaly [Eq. In Eq. ([5]), Aq plays a similar role. 

Since these equations are expressed in terms of the mass 
of the particle consisting of the medium and its number 
density, one can extend them to genuine hadronic matter 
as 



M had 
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M had = £ p . m . A * G . 
i— hadrons 



(G) 
(7) 



The quantities with subscript i denote hadronic counter- 
parts for the nucleon values appearing in Eqs. (jj) and 
([5]). As a simple model, we use a hadron resonance gas 
including all hadrons for which the quantum numbers 
are known as given in the Particle Data Group [27j • The 
number density pi is now generalized to a hadron gas and 
is calculated as a function of T and p,B as 



d l 
2^2 



p 2 dp 



exp[( N /p 2 +™?-^)/T]±l 



(8) 



where the sign is + for fermions and — for bosons and di 
is the degree of freedom of the i-th hadron. We take into 
account the baryon number conservation and strangeness 
conservation and neglect isospin chemical potential for 
simplicity. The strangeness chemical potential p s is de- 
termined from the neutrality condition piSi — [28| . 

In Eq. (J6j> , the masses of hadrons in the chiral limit 
m are needed. Note that strange quark contribution 
from Eq. ([3]) and its off-diagonal counterpart also have 
to be subtracted. Thus we will work within the flavor 
SU(3) symmetric limit. At present, we cannot know all 
of hadron masses in the m u = rrid = m s = limit, es- 
pecially those of the highly excited states. Therefore for 
the masses in the three-flavor chiral limit, we use differ- 
ent masses only for the Goldstone bosons, ground state 
octet and decuplet baryons, and keep the masses of other 
hadrons the same as their vacuum values. Detailed lat- 
tice studies on hadron masses, as done in Ref. [29j with 
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FIG. 1: (Color online) Comparison of the resonance gas with 
lattice results. Upper: interaction measure (e—3p)/T 4 . Lower 
: its gluonic part Mo. Lattice data are taken from the 
HotQCD collaboration for HISQ action with N T = 6 and 8 
[33l ] and from the Budapest-Wuppertal collaboration for stout 
action with continuum estimation (average of N T = 8 and 10) 
[TtJ . For Mo , stout data is estimated from e — 3p by assuming 
the same ratio of the gluonic part as that of HotQCD. See 
text for details. 



physical quark masses, will be helpful for more accurate 



treatment. Specifically, we first put m° = — 0, and 
m° N = 750 MeV from heavy baryon chiral perturbation 
theory 30] ; these are the most important inputs needed 
for the masses in the chiral limit as the contributions 
to the thermodynamic quantities are dominated by these 
hadrons, especially by the Nambu-Goldstone bosons. For 

m„ 



the vector and axial vector mesons, we assume m° p 
and m® = m ai 



We also assume m\ 
m®, rri 



more, we also put mP f - " ,,) " ■" '"" 



to 







= m 



Further- 



K * 

Q 



m\ = rni = = m N , to^» = mh, = rrin = m 
according to the flavor SU(3) symmetry. 

The second moment of the gluon distribution function 
is set to Aq(8to^) = 0.9 for all the hadrons. Generally 
it can differ among hadrons, but it can be shown that 
Aq differs little from this value at such a high energy 
scale where the parametrization of the gluon distribution 
functions [U are relatively well known. 

The result of Mq from the resonance gas model cal- 
culated with Eqs. is shown in the lower panel of 
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FIG. 2: (Color online) Contour plots for Mo [(a)] and M 2 [(b)] 
obtained from the resonance gas model, Eqs. © and ©. The 
thick solid lines (black) indicate the chemical freeze-out line 
including uncertainty of parameters given in Ref. *23]. The 
thick dashed lines (red) also indicate the freeze-out line but 
from [H. 



Fig. [T] together with e — 3p in the upper panel. We 
compare the resonance gas model with lattice data from 
two different fermion discretization schemes. One is from 
highly improved staggered fermion (HISQ) action, calcu- 
lated by the HotQCD collaboration with temporal ex- 
tent iV r = 6 and 8 [33|. The other is from stout- link 
improved staggered fermion (stout) action calculated by 
Budapest-Wuppertal collaboration [l7j . In the former 
the light quark mass is slightly heavier than physical one 
and the continuum extrapolation is not made while the 
latter corresponds to physical quark masses and gives 
a continuum estimation. The upper panel shows the full 
trace anomaly including both gluonic and fermionic parts 
for a reference. Lattice data obtained from the different 
schemes show reasonable agreement in temperature range 
considered here. Therefore we assume the present HISQ 
data already approximates the continuum result well. As 
already discussed in Ref. [13] , the resonance gas model 
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shows a small discrepancy at T > 150 MeV which could 
be attributed to missing heavier states [13] ■ Since in our 
model Mo is essentially dominated by light hadrons, we 
presume that this discrepancy does not affect the follow- 
ing analyses. In the lower panel, the HISQ data shows 
the gluonic part of the trace anomaly. Since the equation 
of state of the stout action was calculated in a different 
way such that the gluonic part is not separated fl7| . we 
estimate it in the following way. First we assume the ra- 
tio of the gluonic part of the trace anomaly to the total 
one is same as those in the HISQ data. Next, we cal- 
culate the ratio by averaging that of HISQ data for 130 
MeV < T < 170 MeV. The upper bound corresponds 
to upper limit of the pseudocritical temperature, below 
which we do not see clear temperature dependency in 
the ratio. Then we multiply e — 3p in the stout action 
by the resultant ratio factor 0.667. Errors are estimated 
from the upper and lower value of the data points and 
from the average deviation of the ratio factor 0.052 of 
the HISQ data. One sees our model M reproduces the 
lattice data well, and therefore we expect that the model 
gives a good approximation to M2 2 We display Mo and 
M2 as functions T and (13 in Fig. [5] 

In Fig. [21 we also draw the chemical freeze-out lines 
proposed by two groups. One (denoted by "FOP) is 
from a combined fit to statistical model results and has 
been shown to agree with various freeze-out conditions 
[23l | . The temperature is given by 



T(fj, B ) = a - bfi% - cfi% 



(9) 



where a = 0.166 ± 0.002 GeV, b = 0.139 ± 0.016 GeV" 1 , 
and c = 0.053 ± 0.021 GeV -3 . Collision energy depen- 
dence is also given through the chemical potential 



161 ± 4 MeV, a' = 1303 ± 120 MeV 
and b' = 0.286 ± 0.049 GeV" 1 . We draw two lines for 
each freeze-out curves corresponding to upper and lower 
temperatures estimated by uncertainties in the param- 
eters. The main difference between the two freeze-out 
curves appears \lb > 300 MeV, which corresponds to 
y/SNN < 10 GeV. While the one by Cleymans et al., 
FOI, has a strong curvature which leads to lower freeze- 
out temperature in this region, the other by Andronic et 
al., FOII, shows almost constant freeze-out temperature 
up to hb ~ 500 MeV, resul ting in coincidence with the 
QCD chiral transition line 37] 3 . Although the difffcr- 
ence between them seems to partly come from the fact 
that FOI uses 4tt particle yield while FOII uses midrapid- 
ity data only, one may consider the constant temperature 
case (FOII) to be simultaneous chemical freeze-out at the 
hadronization. 

For charmonium production, it is not clear which sce- 
nario is more likely, due to small production rate at lower 
energies. At the top SPS energy where the charmonium 
particle ratio data is available, ^/saT/v = 17.3 GeV, the 
two freeze-out curves coincide. From Fig. [2] one sees 
both Mo and M2 increases as \ib does so at fixed tem- 
perature. This implies larger medium modification of 
charmonia for larger chemical potential. If the freeze-out 
temperature decreases steeper, as in FOI, however, re- 
sultant Mo and M2 do not differ so much. Therefore, we 
expect that the mass shift of charmonium substantially 
differs at lower collision energies between the two possible 
freeze-out scenarios. 



III. MASS SHIFT OF CHARMONIUM 



1 + e^/SNN 



(10) 



with d = 1.308±0.028 GeV and e = 0.273±0.016 GeV -1 . 
It has been shown that this parameterization works well 
for recent STAR data [Hj]. 

The other (denoted by FOII) is a paramctrization of 
results of a statistical model shown in Ref. [24]. The 
freeze-out temperature and chemical potential are given 
as functions of ^snn (in unit of GeV) 



r(V5) = T lim ( 1 - 

Mb(Vs) 



1 



0.7+(ev^-2.9)/1.5 



1 + b'y/s NN 



(11) 



(12) 



A. Second order Stark effect 

First, we calculate the mass shift of J/tp using the sec- 
ond order Stark effect in QCD as done in Ref. [1_3]. Pro- 
vided the wave function of the quarkonium in the mo- 
mentum space ip(k) is normalized as J \ip(k)\ 2 = 1, 
the formula of the mass shift for the IS state is given by 



1 

18 



kdk 2 



k 2 /m c 



dip(k) 



7ir 2 a 2 I a s . < 

( —AE 

18 e \ 7T 



dk 



—AE' 

7T 



(13) 
(14) 



2 One may try to improve the agreement by introducing interac- 
tions in the resonance gas. For example, we can use an excluded 
volume correction to the resonance gas, which incorpolates the 
repulsive interaction among hadrons |35|I . We found, however, 
that x 2 fitting of the excluded volume parameter vo to the lat- 
tice data of Mq and e — 3p gives a consistent result with vq = 0. 



where k — \k\, m c and e are the charm quark mass and 
the binding energy, respectively. The above formula can 



3 Although this line corresponds to the chiral transition, it can 
presumably represent the deconfinement transition line also. 
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FIG. 3: (Color online) Temperature dependent part of the 
electric condensate {^-AE 2 }. Each line stands for the case 
of different chemical potential. 

be also derived from potential non-relativistic QCD (pN- 
RQCD), as shown in the Appendix. The second line is 
obtained for the Coulombic bound state with Bohr radius 
a. These parameters can be determined by a fit to the 
J I tp mass in vacuum and the size of the wave function in 
the Cornell potential model p|. It gives m c — 1704 MeV, 
a = 0.271 fm, and a s = 0.57. In this formula, the mass 
shift is proportional to the change of the electric conden- 
sate (^A£ 2 )r ifIB from its vacuum value. The electric 
condensate as well as the magnetic counterpart can be 
written in terms of M and M 2 as [Hj], for Nt = 3, 

(^AE 2 ) = 2 M (T^ B ) + ^M 2 (T^ B ), 

(15) 

/— A£? 2 \ = ~Af (7> B ) + ^M 2 (T, MS ). 

\ 7T / T,hb " 4 7T 

(16) 

Here the effective coupling constant af* can be chosen 
according to the relevant energy scale to the expectation 
value of the operator. In this case, the formula is based 
on OPE with separation scale e. Thus it is plausible to 
take a° s s = 0.57 obtained from the fit to the bound state. 

Figure [3] shows the electric condensate for fi B = 
0, 100, 200, 300 and 400 MeV as a function of tempera- 
ture. The maximum chemical potential 400 MeV roughly 
corresponds to 40.4 GeV Pb+Pb collisions at SPS, the 
lowest collision energy above J/tp production threshold 
[4~fl | . One sees larger change of the electric condensate for 
high temperature and chemical potential as is expected 
from Fig. [21 Now we are able to estimate the mass shift of 
J /if) by putting the condensate into Eq. (fP3"]) . We com- 
pute the mass shift along the freeze-out lines shown in 
Fig. [2] and show the result as a function of ^/snn in Fig. 3] 
For higher colliding energies than ^snn > 30 GeV, the 
mass shift is independent of the colliding energy owing to 
the fact that the freeze-out temperature varies little and 



FIG. 4: (Color online) Mass shift of J/i/j at freeze-out tem- 
perature and chemical potential from the second order Stark 
effect. The horizontal axis denotes the collision energy which 
is related to the temperature and the chemical potential via 
Eqs. @ and (JTUJ) for the FOI case and Eqs. ITT]) and (TT2]) for 
the FOII case. 

the chemical potential does not change the condensate 
significantly (see Fig. [3]) The amount of the downward 
mass shift is 10-20 MeV, including uncertainty. On the 
other hand, low energy results differ between FOI and 
FOII, as expected from Fig. [2] Along the FOI chemical 
freeze-out line, the mass shift becomes smaller while op- 
posite behavior is seen for FOII. In the FOII case, since 
temperature is almost constant, the larger the chemical 
potential, the bigger the mass shift becomes owing to the 
effects from the chemical potential. In the other case, 
however, decreasing freeze-out temperature cancels the 
effect of the chemical potential. For instance, one can 
see in Fig. 0] that AE 2 at T = 165 MeV and fi B = 
is almost equal to that at T = 145 MeV and [is = 400 
MeV. Therefore, At the lowest SPS energy, the down- 
ward mass shift ranges from 10-60 MeV, depending on 
the choice of the thermal parameters. 

B. QCD sum rules 

While the second order Stark effect provides the down- 
ward mass shift directly in the case of the increasing elec- 
tric condensate, QCD sum rule is expected to be more 
quantitatively reliable according to the larger separation 
scale by going to the deep Euclidean region. Here we give 
an estimation based on the Borel sum rule framework 
used in Ref. [l4j with further improvement as described 
below. 

The two point current correlation function in the vec- 
tor channel after the Borel transformation with the Borel 
mass M 2 is given by [IH liH 

M{M 2 ) = e- v TrA(v)[l + a s {M 2 )a{v) + b(u)<h(T) 

+ c(i/)0 c (T)]. (17) 
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with v = Am 2 /M 2 . Here the system is assumed to be at 
rest with respect to the medium. The Wilson coefficients 
a(v), b{v), and c(v) are listed in Ref. [HJ. The tempera- 
ture dependency is governed by the dimension four gluon 
condensate terms <fib{T) and (f> c (T) given by 



An 2 



9(4m2) 2 

4tt 2 
3(4to2)2 



Go(T), 
G 2 (T), 



(18) 
(19) 



The relations of the positive definite quantities Mq and 
M 2 to the dimension four gluon condensates appearing 
in the OPE side correlation function are given by 

G (T) = - ~M (T) 
a cS 

G 2 (T) = -^«-M 2 T 

7T 

after taking the one-loop expression for the beta func- 
tion. As in the pure gauge case, we use a tempera- 
ture dependent effective coupling constant a* = a qq (T) 
extracted from lattice calculation of the color singlet 
heavy quark free energy by assigning a s {G a pp G aup )T = 
(a s (T)G pp G al/p ) in the spirit of the separation scale in 
the heavy quark system which imposes all the tempeture 
effect on the condensates [IH, H2|. We take the values 
from Nf — 2 results of a qq (r max ) in Fig. 6 of Ref. [il j. 
To account for different critical temperatures between the 
Nf = 2 simulation, T c — 202 MeV, and the reality (see 
Sec. H|, we rescale the coupling constant read off from 
the data by assuming similar temperature dependency 
to that of T c = 170 MeV. This value is considered to be 
the upper bound of the pseudocritical temperature in re- 
ality, owing to measurement based on the strange quark 
number susceptibility T pc = 165(5) (3) MeV in Ref. [H 

Presumably, this choice does not affect the results 
quantitatively since the twist-2 contribution to the 
medium modification is relatively small in the hadronic 
phase [ll|. The T = part of the scalar gluon conden- 
sate Gg ac is fixed to be (0.35GeV) 4 [45| as in our previous 
calculations. 4 

The Borel-transformed correlation function is related 
to the spectral density through the dispersion relation 



/>oo 

M(M 2 ) = dse- s/M2 lmfl{s). 
Jo 



(20) 



We model the right-hand side of the dispersion relation 
with a simple ansatz and call it the phenomenological 
side as usual. In the hadronic medium below T c , pre- 
vious analyses based on the gluon condensate of pure 



has still a large error after fitting to the various experimen- 
tal data, see also I4al for example. This vacuum value, however, 
does not affect the in-medium effect significantly since we are 
looking at relative changes from vacuum. 



gauge theory indicates the broadening is small enough to 
ignore. Provided the continuum part of the model spec- 
tral density Ai cont (M 2 ), the the mass of J/-0 is given 
by 



in 



J/4> 



(M 2 ) 



-Q^jjjTjlMjM 2 ) — M cont (M 2 )] 
M(M 2 ) - M cont {M 2 ) 



(21) 



We use the perturbative expression up to 0(a s ) for 
M cont (M 2 ) as in Ref. QJ]. 

Since the mass of J/tp is a function of the Borel mass 
M 2 which is an unphysical parameter, one has to choose 
the range of M 2 called Borel window by following the 
criteria; 

1. M^ in : Convergence of the OPE by imposing the 
dimension four operator contribution less than 30% 
to the total OPE [13]. 

2. M max : Continuum contribution to the dispersion 
integral is less than 30%. We choose the threshold 
parameter s$ such that extracted J / ip mass is least 
sensitive to M 2 . 



As discussed in Ref. [14|, the values 30% are physically 
reasonable but arbitrary. Due to truncation of the OPE, 
we cannot obtain the completely M 2 independent mass. 
Specifically the mass strongly varies with M 2 at lower M 2 
even inside the Borel window. As this can be regarded 
as a systematic uncertainty due to the truncation, we 
take this effect into account in the mass evaluation by 
averaging the mass over the Borel window and take its 
variance as the error [22j]. Namely, 



m 



dM 2 m(M 2 ) / (Af£« - M»J (22) 



and 



{8m) 2 = 



dM 2 (m{M 2 )-rh) 2 / (M, 



2 

max 



M- 



(23) 

An example taken from T — is shown in Fig. [5J Here 
the Borel window is defined by M 2 G [M 2 '- m , M n 2 lax ] and 
M* n ^ma X (Ml in ,M 2 ). 

We have introduced Mq such that dm(M 2 :^/so = 
oo )/dM 2 = in order to remove the strongly M 2 depen- 
dent part of m(M 2 ) from the evaluation of the average 
(|2"21 and variance (|2"3"|) . This is a reasonable choice as the 
continuum threshold is so determined that it makes the 
Borel curve flattest at M 2 > Mq. 

Figure [5] shows an example of the determination pro- 
cess. We start with ^/sq — oo case which gives Mq shown 
as the dashed line. Then we search for ^/sq such that it 
gives the smallest 5m using Eq. (|2"3")l which takes its min- 
imum when the deviation from the average value is the 
smallest. The resultant average deviation is indicated 
by the band in the figure. Irrespective to the tempera- 
ture, it is found to be approximately 5 MeV and 6 MeV 
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3.3 



„ 3.2 
> 

O 



3.1 



\\ Sq =3.49 GeV — 




m=3095 +5MeV 



JlMf, T=0 



1.5 



2.5 



M 2 [GeV 2 ] 



FIG. 5: (Color online) Borel curves at T = 0. The dashed line 
stands for case of yfso = oo. The solid line denotes the case 
of y^so = 3.49 GeV which gives the flattest curve according 
to Eq. (I23|l . The band indicates the systematic uncertainty 
associated with the flattest Borel curve. 



in the case of J/tp and Xci, respectively. Parameters of 
the theory are fixed to m c (p 2 = — 2m 2 ) = 1.262 GeV 
and a s (8m 2 ) = 0.21 by fitting to the vacuum J/ip and 
Xci masses on the basis of the same criterion of the Borel 
window. This process removes the ambiguity on the arbi- 
trary choice of the criterion mentioned above. Including 
the width is straightforward. As shown in Ref. [lij], in- 
troducing width increases the mass at small M 2 . This 
fact leads to larger sq after minimizing dm. Then we 
will have the mass- width relation similar to those shown 
in Refs. pH ED, El- In most cases M§ > M^ in holds 
in the charmonium sum rules. At high temperature and 
chemical potential, however , we found that the Borel sta- 
bility is lost [13, G3] thus Mq is not well denned. In such 
cases, we can still recover the Borel stability by decreas- 
ing the threshold parameter or by introducing the width. 
When the width must be introduced, we cannot deter- 
mine both mass and width simultaneously but have only 
constraints. In what follows, we restrict ourselves to cases 
in which the Borel stability is established with vanishing 
width. 

Figure [5] displays the results of the mass shift obtained 
from the QCD sum rule analysis as described above. We 
plot the mass shifts corresponding to the two freeze-out 
curves as in the Stark effect results (Fig.|4]). Errors are es- 
timated from the uncertainty in the thermal parameters 
in Eqs. (!§))- (TT2"1) . The systematic errors in the Borel sum 
rules which have been introduced above are not included 
in the plot. We also calculate mass shift of Xci in the 
same way. In Xc, we have observed the loss of the Borel 
stability at lower collision energies than ^snn < 8 GeV 
in FOII, owing to much change of the gluon condensate 
1 1 ij . Comparing the result with Fig. |4l one finds that the 
two methods giv e consistent mass shifts as having been 
found in Ref. [14j j . The mass shift of Xci is approxim ately 
twice as large as that of J/-0, as previously found [13.114). 
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FIG. 6: (Color online) Mass shift of J/tp (upper) and \ c \ 
(lower) obtained with QCD sum rules. Solid and dashed lines 
stand for the mass shift corresponding to the different freeze- 
out curves, as in Fig. [4] Errors are calculated from the un- 
certainty in the thermal parameters. 



We expect similar results for other Xc states. 

Before closing the section, we would like to comment on 
the effect of the scattering term, which was pointed out in 
Ref. [48[ to appear as a pole at the zero mode in the corre- 
lator. In our previous works, it was neglected since such 
a contribution appears in the OPE side to cancel the phe- 
nomenological side put as an ansatz. In the deconfined 
phase, this argument should hold because the physical 
particle absorbing the current is the (anti)charm quark, 
while this is not so in the hadronic phase where charmed 
mesons are the physical particles. In general, these terms 
can be neglected as they contribute at zero energy in the 
spectral density. However, without invoking such argu- 
ments, it can be neglected in the present case on the fol- 
lowing grounds. First, the scattering terms in the OPE 
and the phenomenological side will be proportional to 
e -m c /T aiK j e -m D /T respectively, while the other OPE 
terms in the Borel transformed sum rule will in general 
scale as e - 4m c/ M . Hence as long as T < M 2 /{4,m c ) 
or T < M 2 m,D/('irn 2 .), the scattering terms can be ne- 
glected. Since the the smallest Borel mass relevant in our 
analysis is always larger than 1 GeV, taking m c = 1.26 
GeV, one finds that the scattering terms can be safely 
neglected for T < 200 MeV. Moreover, the open charm 
meson will receive greater medium effect than charmonia, 
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making mo close to m c . Thus the scattering contribution 
from the OPE and the phenomenological side will tend to 
cancel each other. Second, one can remove the scattering 
contribution by making use of the fact that it contributes 
as a constant term to the Borel transformed correlator 
M.{M 2 ) — f dse~ s / M p(s) since the scattering contribu- 
tion to the spectral density takes a form as p scat (s) oc 5(s) 
[48j (eauivalentlv y/s5(y/s) (4lj - [5l"l ]). Therefore, the effect 
of the scattering term does not exist in the derivative of 
M(M 2 ) with respect to M 2 (l/M 2 in practical calcula- 
tions). One may then start from the once differentiated 
sum rule and express the mass in the T — limit as 



w fL^{M(M 2 )-M c ° nt (M 2 )} 
-T^M^lMiM^-M^iM 2 )} 



(24) 



Indeed such a method was advocated in Ref. |52[ but 
criticized in Ref. [53[ because starting with a higher or- 
der derivative make the OPE side more sensitive to un- 
known higher dimensional condensates and Borel stabil- 
ity is lost; it was claimed in Ref. [53| that due to this 
artifact the light vector meson was found to increase in 
the medium in Ref. [52| . This comes from the fact that 
the OPE in the light vector meson has no scale parame- 
ter other than the Borel mass; therefore after the Borel 
transformation the sum rule becomes a polynomial in 
1 /M 2 with the highest power determined by the highest 
dimensional operator calculated in the OPE. In the case 
of heavy quarkonia, however, the presence of heavy quark 
mass does not make the OPE a mere polynomial in 1 /M 2 
so that stability is not lost even after derivatives. 

Figure [7] shows an evidence for the above argument. 
We show the Borel curves obtained from Eq. as 
well as those from the ordinary method, Eq. (|2ip . For 
illustration, we show a case with large medium effect, 
corresponding to ^/snn — 8.7 GeV in FOII, of which 
temperature and chemical potential are T = 156 MeV 
and he — 403 MeV, respectively. For vacuum, we dis- 
play three curves. The thin dotted (black) shows the 
same one as in Fig. [5] for the reference. The thin solid 
(green) denotes that obtained from Eq. (|2~4"|) with the 
same value of y/so. One sees both curves give almost the 
same mass. After minimizing 5m in the differentiated 
sum rule [Eq. (|2~4")) ]. one gets a slightly smaller mass in- 
dicated by the thick solid (red) line. One notes Eq. (j2"4")l 
gives smaller mass at small M 2 , as expected from the 
fact that it becomes sensitive to higher dimensional op- 
erators. The small discrepancy of the mass can be at- 
tributed to the dimension six contribution which slightly 
increases the mass [54). The remaining two curves are 
for the medium. The thick dotted (red) is obtained from 
the ordinary sum rule while the thick dashed (blue) is 
from the derivative with the optimization. If one com- 
pares them with the corresponding results for vacuum, 
one finds that the mass shifts are almost the same be- 
tween the two sum rules. Hence, we conclude the gen- 
eral properties of the in-medium modification of heavy 
quarkonia at low temperature up to near T ~ 160 MeV 
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FIG. 7: (Color online) Borel curves for the mass obtained 
from the ordinary QCD sum rule and from the differentiated 
one, Eq. (|24p . See text for detailed explanation. 



will not change by including the scattering contribution. 



IV. IMPLICATION FOR EXPERIMENTS 

We have discussed mass shifts of J/tp in hot medium, 
which could be produced in heavy ion collisions. Pro- 
duction mechanism of J/rp has not been fully understood 
yet, because of the still unknown elementary production 
process and also the complicated collision processes [55[ . 
In the following we briefly review the collision process 
and specify the situation we will consider in this section. 
Charmonia have been considered to be mostly produced 
by collisions between initial state quarks and gluons in 
nuclei at the initial stage of heavy ion collisions. The for- 
mation time scale r ~ 1/ (2m c ) is supposed to be shorter 
than the thermalization time scale of the medium, which 
is related to flow measurements through hydrodynamic 
model calculations. The produced charmonia will also 
interact with colliding nuclei. The dissociation of J/ip 
by this interaction, called cold nuclear matter effect, is 
estimated by the nuclear absorption cross section, which 
is roughly 1.5 mb at RHIC energies and 4.4 mb at SPS 
energies [56j. In the hot medium, charmonia could melt. 
While Lattice QCD have shown existence of the spec- 
tral peak even at higher temperature 0, [EBB, model 
calculations can explain the lattice data with melting of 
J J ip [571 ] . Even if the bound states can survive, they will 
acquire substantial collisional broadening thr oug h inter- 
acting with quarks and gluons in medium [55ll58l| . There 
could be also recombination of a cc pair inside deconfincd 
medium below dissociation temperature. Finally, (anti- 
)charm quarks hadronize at the phase boundary to form 
charmed mesons, baryons, and hidden charm states. If 
J/tp mass is modified in the medium and decay inside 
the mediun, one may be able to observe it as modifica- 
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tion of the peak in dilepton channel. However, one needs 
a dynamical approach in order to take into account var- 
ious process described above, as well as microscopic in- 
formation like cross sections to estimate the final yield 

BUI, Si- 
Here we consider an alternative possibility of indi- 
rect observation via statistical production [H, ||l|, ■ 
This has been already considered as a part of the contri- 
bution in the transport approach and could give substan- 
tial contribution to the final yields at high centrality [55[ . 
In Ref. [4l| . in- medium effect on the charmonium yields 
was considered as a result of D meson mass modification 
and charm conservation. As emphasized in the litera- 
ture, charm conservation plays an important role in the 
statistical description of charm quarks. However, as for 
the particle ratio between charmonia, the fugacity factor 
cancels and the ratio can be expressed as that of ther- 
mal number densities at given temperature and chemical 
potential. Therefore, we can focus on observables domi- 
nated by the statistical production. 

First we examine N^,/ /jVjM at midrapidity which was 
also investigated in Ref. [25[ , since there was experimen- 
tal data for Pb+Pb collisions at ^snn — 17.3 GeV in 
CERN-SPS. The temperature and chemical potential in 
the two freeze-out curves coincide at this enegy as seen 
from Fig. H We adopt T = 160 MeV and [Ib = 240 
MeV in the following [25[ and discuss possible effects of 
charmonium mass shifts. 

Although the QCD sum rule method cannot assess in- 
medium modification of ip' 5 , it is expected to be strongly 
affected [63[. While applicability of the formula of the 
second order Stark effect in Ref. 38] is questionable for 
the physical ip' , a crude estimation based on the dipole 
nature might be possible. If we assume the mass shift 
scales with the size of the wave function, the mass shift 
of ip' is a factor 4.2 larger than that of J/ip- For the SPS 
data at y/sJjN — 17.3 GeV which we will analyze below, 
Am^pi becomes 63-119 MeV by adopting the result from 
QCD sum rules and fully taking the errors into account. 
Here, we vary ip' mass shift in a broader range than the 
above estimation and calculate the ratio as a function of 
Am,/,/ . The total number of J/ip is obtained by summing 
up decay contributions from ip' and Xcj(J = 0, 1, 2). We 
assume the mass shifts of different % c states to be the 
same as that of Xci and that the branching ratios to 
be the same as their vacuum values. Specifically, the 
branching ratios of ip' and Xcj{J — 0, 1, 2) to J/ip are 
taken to be 0.595, 0.0116, 0.344 and 0.195, respectively 

We plot the result for T — 160 MeV and p, B = 240 



5 In the Borel transformed sum rule, we can obtain the same J/ip 
mass even if we incorpolate the t/>' contribution in the model 
spectral function explicitly Therefore, in- medium effect on 

ip' is incorpolated in the effective threshold parameter -^/sq. This 
shows the strong sensitivity of the sum rule to the lowest pole of 
the spectral function. 
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FIG. 8: (Color online) Particle number ratio at T = 160 
MeV and [ib = 240 MeV as a function of ip' mass shift. The 
band in the lower panel indicates the experimental data taken 
from Pb+Pb collisions at ^/snn = 17.3 GeV measured by 
NA50 collaboration [641 ]. The solid lines stand for the results 
without any mass shift in the all charmonia for comparison. 



MeV together with experimental data taken from Pb+Pb 
collisions at ^snn = 17.3 GeV [64J in the lower panel 
of Fig. [5] One sees that a strong variation of the ratio 
as a function of ip' mass shift. The band for our result 
corresponds to the uncertainty in the mass shift shown in 
Fig. |51 If the ip' downward mass shift is smaller than 100 
MeV, which is a reasonable value according to the crude 
estimation, the expected mass shifts are consistent with 
the experimental data. At present, however, it seems 
difficult to draw a conclusion from this observation; in- 
medium modification of ip' is still controversial. While 
the ground state, J/ip, seems to exist in the vicinity of 
the transition, ip' can be dissociated at finite temperature 
even below transition temperature according to potential 
model calculations [9j . Even if ip' is statistically produced 
with the large mass reduction, it can still be dissolved 
into charmed mesons in hadronic medium thus observed 
yields might be smaller than what was produced at the 
hadronization: the dissociation cross section of J/ip in 
hadronic matter is expected to be smaller than that of 
ip' due to its smaller size [59(. In this case, there must be 
of course enhancement of charmed mesons due to charm 
conservation. The analysis might become complicated if 
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initially produced charmonium could survive the quark- 
gluon plasma phase and the hadron phase. Then there 
will be modification to statistical model prediction as dis- 
cussed in a two component models of charmonium pro- 
duction in a heavy ion collision [59l l65j . However, even 
within the two component model, the thermal produc- 
tion will be dominant at LHC and production ratios can 
reveal vital information at the hadronization point. 

Although experimentally challenging, more promising 
observable could be a ratio of Xc to J/ip. In the up- 
per panel, we display a fraction of J/tp coming from Xci 
and Xc2 to the total number of J/tp for T = 160 MeV. 
This quantity, R Xc was measured by_ HERA-B Collabo- 
ration in proton- nucleus collisions [66j . Our result shows 
10-20% increase of this quantity for the statistical pro- 
duction with mass reductions, almost independent of the 
uncertain tp' mass shift. 



V. SUMMARY 

In this paper, we investigate the mass shift of charmo- 
nia induced by change of gluon condensates in hadronic 
medium by making use of both perturbative (QCD sec- 
ond order Stark effect) and non-perturbative (QCD sum 
rule) approaches. The inputs for the medium effect is the 
gluon condensates calculated by a resonance gas model 
that well reproduces the thermodynamic quantities cal- 
culated on the lattice QCD result in the temperature re- 
gion considered here. Extending it to the finite baryonic 
chemical potential, we found that the change of the gluon 
condensates becomes larger at large chemical potential, 
which could result in a larger medium effect on charmo- 
nium production in heavy ion collisions at lower colliding 
energies, depending on hadronization temperature. We 
found that both the perturbative and non-perturbative 
estimations give almost the same mass shift. After elabo- 
rating on the sum rule analysis, we estimated systematic 
error on the obtained mass shift and found that most 
of the uncertainty comes mainly from thermal parame- 



ters. While the mass shift is almost constant along the 
chemical freeze-out line of which freeze-out temperature 
decreases as chemical potential increases (FOI), it can 
exhibit stronger downward shift at lower colliding ener- 
gies if the charmonia are produced at hadronization and 
simultaneously frozen out (FOII). 

We consider experimental implications of this observa- 
tion in the context of statistical hadronization picture in 
which hadronization and the freeze-out of charmonium 
occur simultaneously. We found the data in Pb+Pb col- 
lisions at 17^4GeV are consistent with mass reduction of 
J/tp and Xc although the effect on ip' should be clar- 
ified before final conclusion is made. We pointed out 
that 10-20% enhancement of the production ratio be- 
tween Xc and J I tp could be a signature of the downward 
mass shift, which indicate a precursor of the deconfinc- 
ment phenomenon. 
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Appendix: Derivation of the second order Stark 
effect formula from Potential NRQCD 

In this appendix, we show that the second order Stark 
effect, Eq. (fT3"|) , can be derived from potential non- 
relativistic QCD (pNRQCD) which provides a systematic 
perturbative approach to the OPE. 



J 



Here the effective lagrangian in the static limit is given by [67|, |68[ , 



n f r f 



id + C> 



VkTV-fOV • gES + S^f-gEO} + -y-Tr{OV • gEO + O f Or • gE} 



a Vs 



S + O 



1 a 



Vo 



2N n r 



O 



(A.l) 



The fields S = Sl c /\fWl and O = O a T a /^/Tp are normalized static quark-antiquark singlet and octet fields, re- 
spectively, E is the chromoelectric field, and Cf — (AT? — 1)/(2N C ) — 4/3. The trace is over the color indices. The 
matching coefficients at leading order are, ay s = ct s ,av = ct s , Va — KVs = 1- The leading order correction to the 
singlet potential at finite temperature as given in Eq. (64) of Ref. |69| is, 



[SVs(r) ]n = -„ 2 g^ 



dte~ ltAV 



E a (t)cP(t,0) ab E b (0) 



(A.2) 
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where 



AV =- 

r 



2N r 



CfO-Vs 



2r 



(A.3) 



and is the potential difference between the singlet ground state and the octet state excited by the color electric field. 
d is the number of the dimension in the regularization of the momentum integral. In Ref. [69| , the electric propagator 
is calculated in thermal perturbation in various limits. 

To obtain the formula for the 2nd order Stark effect Eq. (|A.2|1 . we have to take the following steps. 



1. First, to extract the contribution from the lowest dimensional operator, we take, 



E a (t)4>{t,0) ab E b {Q) 



E a {0)E a {0) 



(A.4) 



Moreover, since we are interested in temperatures near T c , instead of using thermal perturbation to calculate the 
temperature dependent part of the electric condensate, we use the nonperturbative value extracted from lattice 
QCD. This approximation is valid as long as the scales in the matrix element is smaller than the separation 
scale. 

2. We take the matrix element of Eq. (|A.2[) for the ground state charmonium, and calculate it using the relative 
momentum between the cc quarks. 



• I i f d 3 p r- ..-J/ 



SE <» = -^W—l J (2#^^I^WI^^(0)^(0))t (A.5) 

3. We assume that the energies for the intermediate octet charmonium sate and the initial ground state can be 
written as follows, 



Ej/^ = 2m c - e 

E = 2m c +p 2 /m Cl (A.6) 

where e is the binding energy for the J ftp and Eo represents the energy of the octet continuum state. Putting 
these energies into Eq. (|A.5[) and provided Tp = 1/2, d — 4 and N c = 3, one obtains Eq. (|T5|) . The contribution 
from the repulsive coulomb potential was neglected in the continuum octet energy as it vanishes in the large N c 
limit taken in the Peskin formalism [§8j]. 



The leading order OPE term taken here is quite sim- 
ilar to the limit taken in Refs. Iz3|. The difference 
there was the assumption of Lorentz invariance of the 
vacuum which is broken at finite temperature; therefore 
the mass shift was proportional to the gluon condensate. 



As we have identified the approximations taken in the 
derivation of the 2nd order Stark effect, it would be use- 
ful to improve the formula by taking into account the 
renormalization group improved potentials and the 1 /N c 
corrections [681 ] . 
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